/** Fichier d'implémentation de la classe Solveur.
 *
 * Le fichier contient la classe Solveur, qui peut résoudre exactement des
 * équations de degré inférieur ou égal à 4.
 *
 * @file Solveur.cpp
 * @author Johan "Solo" Luisier
 *
 * @date 2009
 */
#include "Solveur.hpp"

#include <algorithm>
#include <cmath>

#include <iostream>

using namespace std;

/** Constructeur par défaut.
 *
 * Constucteur, qui initialise tous les coefficients de
 * l'équation. L'ordre et les valeurs par défaut sont choisis de telle
 * sorte que le degré de l'équation est déterminé par le nombre
 * d'arguments donnés.
 *
 * @param p0 terme constant.
 * @param p1 terme linéraire.
 * @param p2 terme quadratique.
 * @param p3 terme cubique.
 * @param p4 terme quartique.
 */
Solveur::Solveur( const double& p0, const double& p1, const double& p2,
                  const double& p3, const double& p4 )
    : a0( p0 ), a1( p1 ), a2( p2 ), a3( p3 ), a4( p4 )
{}

/** Destructeur, libère la mémoire utilisée.
 *
 * Le destructeur n'a rien de spécial à faire.
 */
Solveur::~Solveur()
{}

/** Méthode qui permet de dàterminer la valeur des coefficients.
 *
 * Méthode permettant de changer la valeur des coefficients. A nouveau
 * l'ordre et les valeurs par défaut permettent d'obtenir un degré
 * quelconque.
 *
 * @param p0 terme constant.
 * @param p1 terme linéraire.
 * @param p2 terme quadratique.
 * @param p3 terme cubique.
 * @param p4 terme quartique.
 */
void Solveur::setParams( const double& p0, const double& p1,
                         const double& p2, const double& p3,
                         const double& p4 )
{
    a0 = p0;
    a1 = p1;
    a2 = p2;
    a3 = p3;
    a4 = p4;
}

/** Méthode qui résoud l'équation déterminée par les coefficients entrés.
 *
 * Interface principale, effectue la résolution en appelant une
 * méthode spécifique en fonction du degré de l'équation.
 *
 * Les solutions obtenues sont triées par ordre croissant avant d'être
 * renvoyées.
 *
 * @param solution conteneur dans lequel seront stockés les
 * éventuelles solutions.
 */
void Solveur::resoud( vector< double >& solution )
{
    solution.clear();
    if ( a4 == 0. && a3 == 0. && a2 == 0. )
        resoudDegre1( solution );
    if ( a4 == 0. && a3 == 0. )
        resoudDegre2( solution );
    else if ( a4 == 0. )
        resoudDegre3( solution );
    else
        resoudDegre4( solution );

    sort( solution.begin(), solution.end() );
}

/** Méthode de résolution des équations de degré 1.
 *
 * Méthode qui résoud les équation de degré 1. Ne sert à rien, mais
 * comme ça l'outil est complet.
 *
 * @param solution conteneur dans lequel sera stockée l'éventuelle
 * solution.
 */
void Solveur::resoudDegre1( vector< double >& solution )
{
    if ( a1 == 0. )
        return;

    solution.push_back( - a0 / a1 );
}

/** Méthode de résolution des équations de degré 2.
 *
 * Méthode qui résoud les équations de degré 2 par la méthode
 * habituelle du discriminant. Les nombres dont la valeur absolue est
 * inférieure à 1e-7 sont arrondis à zéro.
 *
 * @param solution conteneur dans lequel seront stockées les
 * éventuelles solutions.
 */
void Solveur::resoudDegre2( vector< double >& solution )
{
    double Delta( a1*a1 - 4 * a0 * a2 );

    if ( Delta < -1e-7 )
        return;
    else if ( abs( Delta ) < 1e-7 )
        solution.push_back( - a1 / ( 2 * a2 ) );
    else
    {
        solution.push_back( ( - a1 + sqrt( Delta ) ) /( 2 * a2 ) );
        solution.push_back( ( - a1 - sqrt( Delta ) ) /( 2 * a2 ) );
    }
}

/** Méthode de résolution des équations de degré 3.
 *
 * Méthode qui résoud les équations de degré 3 par la méthode de
 * Cardan. Les nombres dont la valeur absolue est inférieure à 1e-7
 * sont arrondis à zéro. La méthode repose sur le calcul d'un
 * déterminant.
 *
 * L'équation est en premier lieu normalisée pour que le coefficient
 * dominant soit 1. Il ne peut y avoir de problème de division par
 * zéro, car on serait en présence d'une équation du 2ème degré.
 *
 * @param solution conteneur dans lequel seront stockées les
 * éventuelles solutions.
 */
void Solveur::resoudDegre3( vector< double >& solution )
{
    a0 /= a3;
    a1 /= a3;
    a2 /= a3;
    a3 /= a3;

    double q( ( a2*a2 - 3 * a1 ) / 9. );
    double r( ( a2 * (a2*a2 - 4.5 * a1) + 13.5 * a0 ) / 27. );
    double D( q*q*q - r*r );

    if ( D > 0. || abs( D ) < 1e-7 )
    {
        double Dp( r / sqrt( q*q*q ) );
        double tau( acos( Dp ) / 3. );
        double s( -2 * sqrt( q ) );
        double PI( acos( -1. ) );
        for (unsigned int i(0); i < 3; i++)
            solution.push_back( s * cos( tau + 2 * PI * i / 3. ) - a2 / 3. );
    }
    else
    {
        double s( pow( sqrt(-D) + abs( r ), 1./3.) );
        if ( r  < -1e-7 )
            solution.push_back( s + q / s - a2 / 3. );
        else
            solution.push_back( - s - q / s - a2 / 3. );
    }
}

/** Méthode de résolution des équations de degré 4.
 *
 * Méthode qui résoud les équations de degré 4 par la méthode de
 * Ferrari.
 *
 * L'équation est en premier lieu normalisée pour que le coefficient
 * dominant soit 1. Il ne peut y avoir de problème de division par
 * zéro, car on serait en présence d'une équation du 3ème degré.
 *
 * Dans le cas où le terme constant est nul, le cas est ramené à une
 * équation de degré 3, modulo le fait que 0 est une solution de
 * l'équation de degré 4.
 *
 * On utilise également le fait que dans certains cas l'équation peut
 * se ramener à une équation bicarrée, que l'on résoud alors en
 * faisant appel à resoudDegre2.
 *
 * @param solution conteneur dans lequel seront stockées les
 * éventuelles solutions.
 */
void Solveur::resoudDegre4( vector< double >& solution )
{
    a0 /= a4;
    a1 /= a4;
    a2 /= a4;
    a3 /= a4;
    a4 /= a4;

    if ( abs( a0 ) < 1e-7 )
    {
        solution.push_back( 0. );
        a0 = a1;
        a1 = a2;
        a2 = a3;
        a3 = a4;
        a4 = 0.;
        resoudDegre3( solution );
    }
    else
    {
        double alpha( -3. * a3 * a3 / 8. + a2 );
        double beta( pow( a3, 3 ) / 8. - a3 * a2 / 2. + a1 );
        double gamma( - 3. * pow( a3, 4 ) / 256. + a2 * pow( a3, 2 ) / 16.
                      - a3 * a1 / 4. + a0 );

        if ( abs( beta ) < 1e-7 )
        {
            Solveur biCarre( gamma, alpha, 1.);

            vector< double > secondaire;

            biCarre.resoud( secondaire );

            vector< double >::const_iterator solIt,
                    solBegin( secondaire.begin() ), solEnd( secondaire.end() );

            for (solIt = solBegin; solIt != solEnd; solIt++)
            {
                if ( abs( *solIt ) < 1e-7 )
                    solution.push_back( 0. );
                else if ( *solIt > 1e-7 )
                {
                    solution.push_back( sqrt( *solIt ) );
                    solution.push_back( - sqrt( *solIt ) );
                }
            }

            /*
        double racine1( trouveRacine( pow( alpha, 2 ) - 4. * gamma ) );

        if ( racine1 < 0. )
        return;
        double racine2( trouveRacine( ( - alpha + racine1 ) / 2. ) );
        if ( racine2 > 0. )
        {
        solution.push_back( - a3 / 4. + racine2 );
        solution.push_back( - a3 / 4. - racine2 );
        }
        else if ( racine2 == 0. )
        solution.push_back( - a3 / 4. );
      */
            return;
        }
        else if ( abs( gamma ) < 1e-7 )
        {
            solution.push_back( 0. );

            Solveur resolvante( beta, alpha, 0., 1. );

            resolvante.resoud( solution );

            return;
        }
        else
        {
            double p( - 3. * a3*a3 / 8. + a2  );
            double q( a3*a3*a3 / 8. - a3 * a2 / 2. + a1 );
            double r( - 3. / 256. * a3*a3*a3*a3 + a3*a3 * a2 / 16. - a3 * a1 / 4.
                      + a0  );
            double b3( 8. );
            double b2( - 4. * p  );
            double b1( - 8. * r );
            double b0( 4. * r * p  - q * q  );

            vector< double > secondaire;

            Solveur aide( b0, b1, b2, b3 );

            aide.resoudDegre3( secondaire );

            if ( secondaire.empty() )
                // pas de solution
                return;

            double y0( secondaire.front() );

            double tmp( - p + 2. * y0 ), A0, B0;

            if ( tmp > 1e-7 )
            {
                A0 = sqrt( tmp );
                B0 = - q / ( 2 * A0 ) ;
            }
            else
                return;

            Solveur eqAux1( y0 - B0, -A0, 1. );
            Solveur eqAux2( y0 + B0,  A0, 1. );

            vector< double > zi;

            eqAux1.resoudDegre2( zi );
            eqAux2.resoudDegre2( zi );

            vector< double >::const_iterator Begin( zi.begin() ), End( zi.end() ), It;

            for (It = Begin; It != End; It++)
                solution.push_back( *It - a3 / 4. );
        }
    }
}

/** Méthode qui traite le calcul de la racine carrée d'un nombre.
 *
 * Méthode qui renvoie la racine carrée d'un nombre, en tenant compte
 * des erreurs d'arrondi : tout nombre inférieur en valeur absolue à
 * 1e-7 est arrondi à zéro.
 *
 * @param carre nombre dont on veut extraire la racine carrée.
 * @return la racine carrée du nombre, ou -1 s'il s'agit d'un nombre
 * négatif.
 */
double Solveur::trouveRacine( const double& carre ) const
{
    if ( carre < -1e-14 )
        return -1.;
    else if ( abs( carre ) < 1e-14 )
        return 0.;
    else
        return sqrt( carre );
}
